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Abstract. Measurements of the lateral distribution of 
Cherenkov photons with the wide-angle atmospheric 
Cherenkov light detector array AIROBICC and of the 
charged particle lateral distribution with the scintillator 
matrix of the HEGRA air-shower detector complex in air 
showers are reported. They are used in conjunction to de- 
termine the energy spectrum and coarse chemical compo- 
sition of charged cosmic rays in the energy interval from 
0.3 PeV to 10 PeV. With the atmospheric shower-front 
sampling technique these detectors measure the electro- 
magnetic component of an extensive air shower via the lat- 
eral density distribution of the shower particles and of the 
Cherenkov photons. The data are compared with events 
generated with the CORSIKA program package with the 
QGSJET hadronic-event generator. Consistency checks 
performed with primary energy-reconstruction methods 
based on different shower observables indicate satisfac- 
tory agreement between these extensive air shower simu- 
lations and the experimental data. This permits to derive 
results concerning the energy spectrum and composition 
of charged cosmic rays. 

The energy spectrum features a so called "knee" 
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at an energy of £;knco=3.98to|^(stat) ± 0.53(syst) PeV. 
Power law fits to the differential energy spectrum 
yield indices of -2.72to;o3(stat) ± 0.07(syst) below and 
-3.22l°;^^(stat) ± 0.08(syst) above the knee. 
The best-fit elongation rate for the whole energy range is 
determined to 78.3 ± 1.0 (stat) ± 6.2 (syst) g/cm^. At the 
highest energies it seems to decrease slightly. The best-fit 
fraction of light nuclei decreases from 37^21% (combined 
statistical and systematic) to 8lg^% (combined statisti- 
cal and systematic) in the energy range discussed here. 
A detailed study of the systematic errors reveals that a 
non-changing composition cannot be excluded. 

Key words: cosmic rays 



1. Introduction 

The origin of extra solar cosmic rays (CR) is one 
of the important unresolved astrophysical questions. 
Galactic shell type supernova remnants (SNR) have 
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been proposed as plausible acceleration sites for cos- 
mic rays up to energies of several PeV (Volk 1997) 
and - for very massive SN progenitors - to even higher 
energies (Biermann et al. 1995). Recently direct experi- 
mental evidence for electron acceleration in these ob- 
jects has been found in the X-ray (Koyama et al. 1995; 
Koyama et al. 1997; Allen et al. 1997) and TeV 7-ray 
range (Tanimori et al. 1998). Somewhat surprisingly, sim- 
ilar searches for evidence of hadron acceleration have 
only yielded upper limits on the expected 7-ray emis- 
sion from the interaction of the hadrons with inter- 
stellar matter up to now (Prosch et al. 1996; Hefi 1997; 
Buckley et al. 1998). 

An indirect approach to distinguish between different 
theoretical models aiming to describe the acceleration 
of charged cosmic rays (CR) is to measure the energy 
spectrum and composition of CR and compare the re- 
sults with model predictions. Here the energy regime 
around the so called "knee" between 1 and 10 PeV is espe- 
cially interesting (Watson 1997). In this energy range the 
all-particle CR energy spectral slope - that is constant 
within measurement errors for lower energies - suddenly 
increases. The riddle of the origin of the knee and of the 
cosmic radiation with energies exceeding it, is not yet fi- 
nally resolved. The following general solutions have been 
discussed: 

1. The change in index is due to some propagation ef- 
fect in an "original" cosmic-ray population that displays 
an unbroken power law from low energies up to energies 
above the knee. The most popular idea is that the en- 
ergy dependence of the diffusion constant of cosmic rays in 
the Galaxy could change in the knee region (Peters 1961; 
Ptuskin et al. 1993). Because of the dependence of the dif- 
fusion constant on the nuclear charge Z a modest decrease 
in the fraction of "light" elements (hydrogen and helium) 
would be expected. In the simplified chemical model we 
use below (heavy elements modelled by 65 % oxygen and 
35 % iron, light elements by 40 % hydrogen and 60 % he- 
lium) the fraction of light elements would be expected to 
decrease from an assumed value of 60 % below the knee to 
43 % above the knee. In such a model (barring a special 
cancellation of effects) the knee is expected to be a rela- 
tively smooth feature, extending over about a decade in 
energy. A principal problem with this approach is that no 
plausible Galactic source of cosmic rays has been identified 
which is quantitatively capable of producing the "original" 
cosmic-ray population. 

2. The knee signals in some way the maximum energy for 
the sources responsible for low energy cosmic rays. The 
cosmic rays at higher energies could be "re-accelerated" 
low energy cosmic-rays, e.g. at the shock front of a Galac- 
tic wind (Jokipii & Morfill 1987) or an ensemble of shock 
fronts in clusters of massive stars (Bykov et al. 1997). In 
this case a phenomenology similar to the diffusion model 
in the previous paragraph would be expected. Alterna- 
tively, above the knee a completely new population of 



cosmic rays dominates. In this case typically dramatic 
changes in chemical composition are expected, e.g. to pure 
hydrogen in the extragalactical model of Protheroe (1992) 
and nearly pure heavy elements (fraction of light elements 
< 0.3 far above the knee) in a model with special SNRs 
by Stanev et al. (1993). The special properties of these 
new sources could in principle allow to understand a knee 
relatively "sharp" in energy. 

To definitely discriminate between a composition changing 
as expected in models with an energy dependent diffusion 
constant (discussed above under 1.) and an unchanging 
composition, it is necessary to achieve an error of smaller 
than ± 10 % in the experimental determination of the frac- 
tion of light elements in the total cosmic radiation above 
the knee. 

While the cosmic-ray composition and energy spectrum 
are well known from direct balloon and space-borne obser- 
vations up to energies of about 100 TeV, no general agree- 
ment has been reached at higher energies (Watson 1997). 
The results obtained for CR around the "knee" suffer se- 
riously from the fact that due to the low flux of CR above 
1 PeV, only large ground based installations observing the 
extensive air showers (EAS) induced by cosmic rays in the 
atmosphere can provide experimental data. However the 
sensitivity of EAS observables to the mass of the primary 
CR is weak. The analyses are rendered even more diffi- 
cult due to theoretical uncertainties concerning the high 
energy interactions in the atmosphere (Knapp et al. 1996; 
Gaisser 1997). 

Here we present an analysis of EAS between 300 TeV and 
10 PeV which restricts to observables related to the elec- 
tromagnetic shower component. In the following sections 
the experimental setup (section 2), the Monte-Carlo sim- 
ulations (section 3), the event reconstruction (section 4) 
and analysis methods (section 5) are described. Sections 6 
presents the results concerning the CR energy spectrum 
and the coarse mass composition. A more detailed study 
of the systematic errors and a discussion of methods to 
analyse the composition without relying on the absolute 
penetration depth are discussed in section 7. The paper 
ends with conclusions in sections. 

2. The experimental Setup 

The air-shower detector complex HEGRA covers an area 
of ISO- 180 m2 at a height of 2200 m 
a.s.l. (790g/cm2) (Lindner et al. 1997; Karle et al. 1995; 
Krawczynski ct al. 1996 ; Rhode et al. 1996). In the 
present analysis only data of the scintillator array and 
part of the AIROBICC array were used. The former con- 
sists of 243 huts with plastic scintillators of an area of 0.96 
m^, covered with 5 mm of lead on a grid with 15 m spac- 
ing (with a denser part with 10 m spacing in the centre of 
the array). The part of the latter used in this analysis is 
formed by 49 open photo-multipliers fitted with Winston 
cones, restricting the viewed solid angle to 0.835 sr and 
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measuring the air Cherenkov light of EAS on a grid with 
30 m spacing. The gain nonlinearity of all components in 
the Cherenkov-light measurement was carefully checked, 
both with a LED light source with variable light inten- 
sity and a direct charge source. While the used Cherenkov 
light photomultiplier tubes were found to be linear, the 
used amplifier showed an antilinearity (gain rises with in- 
put amplitude) which was corrected in the data analysis. 
Above 10 PeV the amplifier begin to show signs of satu- 
ration and therefore no data above 10 PeV are included 
in the present analysis. The trigger conditions used for 
the data analysed here demand a signal from at least 14 
scintillator or 6 AIROBICC stations within 150 ns. This 
corresponds to an energy threshold for primary protons 
and iron nuclei of 25 and 80 TeV respectively. 



3. Monte Carlo Simulations 

EAS events were simulated using the CORSIKA code 
in its versions 5.20 with the QGSJET/GHEISHA op- 
tions (Knapp & Heck 1993; Heck et al 1998)). This gen- 
erator is based on the quark-gluon string model 
(QGS) with an allowance for semihard processes (JET) 
(Kalmykov et al. 1997). Complex nuclei were treated with 
the "complete fragmentation" ansatz. The energy cutoff 
for particles of the electromagnetic cascade was set to 
3 MeV. Proton, helium, oxygen and iron induced show- 
ers were produced at zenith angles of 0°,6°, 12° and 18° 
at discrete energies between 300 TeV and 10 PeV (4400 
independent showers of which 1000 are above 1 PeV) as 
well as an independent sample with a continuous energy 
distribution at zenith angles of 6° and 12° between 50 
TeV and 1 PeV (4330 independent showers) and follow- 
ing a power law of E^^ between 2.5 PeV and 6.5 PeV 
(240 showers). The events continuously distributed in en- 
ergy were spectrally weighted and used in the fits to infer 
the chemical composition (section 5.2), while the samples 
with discrete energies were employed to develop the re- 
constrTiction methods and to correct the results obtained 
with biased estimators of the primary energy (section 5.1). 
Note that the simulation of an EAS induced by a 1 PeV 
primary proton including the Cherenkov light production 
requires about 3 h CPU time on a 300 MHz Pentium-II PC 
(during the same time the HEGRA experiment registers 
more than 350 showers with energies larger than 1 PeV). 
The detector performance was modelled with two in- 
dependent detector simulations: a full detector simu- 
lation (Martinez et al. 1995), and an empirical simula- 
tion using measured response functions ( Haustein 1996; 
Horns 1997). Each independent generated EAS was used 
20 times with core positions inside and outside the 
HEGRA area to take into account the detector related 
fluctuations of observables and to check the event selec- 
tion criteria. With the standard cuts described below, each 
shower was used on the average two times in the discrete 



and once in the continuous MC sample. 
Special care has been taken to simulate the density pro- 
file and absorption features of the atmosphere above 
LaPalma correctly. Weather balloon measurements as 
well as comparisons between TeV photon data registered 
by the HEGRA system of imaging air Cherenkov tele- 
scopes and simulations were employed for this purpose 
(Bernlohr 1999;Konopelko et al. 1999). The shower devel- 
opment and light emission were modelled with the U.S. 
standard atmosphere, and the light propagation was then 
simulated with a special program ELBA (Haustein 1996), 
assuming a tropical maritime atmosphere for the summer. 
This atmosphere is a good approximation for the condi- 
tions at Tenerife, an island neighbouring the experimental 
site (Bernlohr 1999). 

4. Event Reconstruction and Data Selection 

The core position of an EAS is reconstructed indepen- 
dently from the data of the scintillator matrix and from 
AIROBICC where the latter data allow to tag core po- 
sitions beyond the HEGRA boundary. If the core posi- 
tion lies inside the area covered with detector elements 
the scintillator derived core coordinates have a resolution 
of (7(corc) = 2(5) m for protons (iron) at energies above 
300 TeV (a little more accurate compared to AIROBICC 
mainly due to the smaller grid distances of the scintilla- 
tor huts). The direction of the primary particle is recon- 




20 40 60 80 100 120 140 

radial distance (m) 



Fig. 1. The Cherenkov-light (pc, open circles, in 10^ x 
photons/m^) and the actually measured (under lead coverage, 
see text) charged-particlc density (pe, full circles, in 1/m^) as 
determined in a single observed shower. Each open (closed) 
circle corresponds to the light (charged particle) density de- 
termined by one AIROBICC (scintillator) station. The NKG 
function has been fitted to the charged-particlc distribution 
and an exponential function to the Cherenkov light distribu- 
tion in a 7.5 - 90 m radial interval (full lines). The energy of this 
event was 2.1 (2.8) PeV if induced by a proton (iron) nucleus. 
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structed nearly independently from the scintillator and 
AIROBICC arrival time data (where the scintillator de- 
rived core position is used here). 

The particle density measured by the scintillator array is 
fitted by the NKG formula (Greisen 1956) with a Moliere 
radius of 106 m, yielding the shower size Ns and an age 
value. Ng is a factor «1.6 larger than Ng (the "true" 
shower size at detector level, denoting the number of 
charged particles above a kinetic energy of 3MeV) due 
to coverage of the detector huts with a lead layer and the 
fact that the NKG function does not correctly parametrise 
the electromagnetic part of hadronic showers. 
The dependence of the Cherenkov photon density pc, as 
measured by AIROBICC, on the distance r from the 
shower axis can be well described by an exponential in 
the region 20m <r< 100m (Patterson & Hillas 1983): 

/9c(r) = &■ ex.p{i ■ slope) . (1) 

The parameter slope (in units of [1/m]) is the most im- 
portant one in our analysis methods. As an illustration 
Fig.l shows the lateral charged-particle and Cherenkov- 
light distributions for a single shower. 
The amplitude calibration of the scintillator array is done 
for samples of 50000 events by comparing the ADC spec- 
tra of the individual huts - which display a single peak 
essentially corresponding to the energy deposited by min- 
imum ionising electrons and muons - with the result of 
MC simulations for identical conditions. The absolute 
amount of the air Cherenkov light registered by AIRO- 
BICC was calibrated by comparing the energy inferred 
from the lateral Cherenkov light density in the spectral 
range from 300 nm to 500 nm registered at a shower 
core distance of 90 m (referred to as Lgo in the following) 
and the energy derived from Ns and slope in the inter- 
val 2.5 < logio(E(Ns,5/ope)/TeV) < 2.75 (refer to section 
5.1 below for energy reconstruction methods). The abso- 
lute Chcrcnkov-light calibration thus depends on the CR 
mass composition, because we do not apply a primary- 
mass independent energy reconstruction here. We used 
the low-energy composition at 100 TeV as specified by 
Wiebel-Sooth et al.(1998) (60 % light elements, see be- 
low section 5.2 for details) for this calibration. If a pure 
proton (iron) composition is assumed the energy recon- 
structed from Cherenkov light alone is shifted by 3 (13)% 
to higher (lower) energies. 

To select well-reconstructed events the EAS core positions 
and directions as reconstructed with AIROBICC and the 
scintillator matrix are demanded to be consistent. Addi- 
tional cuts ensure the quality of the directional as well as 
the fits to the lateral particle and Cherenkov light den- 
sity distributions. Events with the true shower-core posi- 
tion within the HEGRA array boundaries for the detector 
components used in this analysis (distance to edge of ar- 
ray > 10 111) and a zenith angle below 15° are used for 
the further analysis. The efficiency to select EAS events 
with true core positions in the regarded 160- 160 m^ area 



is about 98% for primary energies above 300 TeV (inde- 
pendent of the primary mass). The contamination of the 
sample with EAS, where the true shower cores lay be- 
yond the HEGRA boundary but which were erroneously 
reconstructed to fulfil the cuts is less than 1% from our 
simulations. 

Nights with perfect weather conditions are selected by 
data of the Carlsberg Meridian Cycle (Argyle 1997) and 
by comparing the Cherenkov light measurements with 
data from the scintillators for samples of 50000 events (ac- 
cumulated in about one hour with the used setup and trig- 
ger conditions). The data set solely contains nights with- 
out any technical problems of the used detector stations. 
In total it comprises (dead-time corrected) an on-time of 
208 h. This corresponds to about 150000 events after all 
cuts with an energy above about 100 TeV and a zenith 
angle below 15°. 

5. The Analysis Methods 

Monte Carlo simulations revealed that the distance d,nax 
between detector and shower maximum (defined as the 
point in the shower development with the maximal num- 
ber of charged particles) can be reconstructed indepen- 
dently of the primary mass with the shape parameter slope 
of the lateral Cherenkov light density distribution: 

di„ax = [680 + slope- 20880 m] g/cm^. (2) 

From this relation the distance to the shower maximum is 
determined with a resolution (i.e. root mean square (RMS) 
of a (dmax(true)-dmaa; (reconstructed) distribution) rang- 
ing from 40g/cm2 at 300 TeV to 20g/cm2 at lOPeV (in- 
cluding all detector effects but no systematic error in the 
mean d^ax)- The most important technical improvement 
in the data presented here to previous experiments is that 
these values are distinctly smaller than the width of nat- 
ural shower fluctuations of proton induced showers in the 
atmosphere (see below Table 3). This makes the shape of 
the penetration-depth distribution a sensitive parameter 
for the chemical c;omposition (see section 7). Simple geo- 
metrical relations permit to infer ATmax, the depth of the 
shower maximum in the atmosphere, from dmax- Relation 
2 is only weakly energy dependent (Lindner 1998a). This 
dependence is neglected here. 

5.1. Energy Reconstruction 

Methods have been developed to reconstruct the pri- 
mary CR energy from the scintillator and AIROBICC 
data independently of the primary mass with an accu- 
racy better than 35% (Lindner 1998a; Cortina et al. 1997; 
Cortina et al. 1998a). However, these methods lead to a 
relatively strong correlation between reconstructed energy 
and Amax (showers with a maxiniuiii position that fluctu- 
ated to smaller values compared to the mean Xmax, are 
reconstructed with higher energies). In order to infer the 
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chemical composition, a careful modelation of the response 
function between the variables Ng (or Lgo) and slope on 
the one hand and energy and composition (or penetra- 
tion depth) on the other hand ( e.g. via two-dimensional 
regularised unfolding) is then necessary. Such procedures 
have been employed in some analyses of HEGRA data 
(Wiebel-Sooth 1998; Kornmayer 1999). Two reasons lead 
us to prefer to circumvent the mentioned problem with 
the use of two simpler energy estimators here. 
One reason is that the methods described below are based 
on physically transparent properties of air-showers in- 
ferred from the Monte-Carlo simulations. Whether these 
properties really hold, is tested to some degree using differ- 
ent energy estimators with different biases and comparing 
the obtained results. These consistency checks are an im- 
portant advantage over more refined and complete meth- 
ods when it is doubtful how well the Monte-Carlo simula- 
tion describes the data. The other reason is that the Monte 
Carlo statistics at the highest energy is still rather lim- 
ited and mean shower properties are inferred with higher 
certainty than a complete response matrix. The mass in- 
dependent energy reconstruction methods will be applied 
to the data in a forthcoming publication together with a 
discussion of the influence of different EAS simulations. 
The energy estimators used in this paper and described 
below are based on Ng and slope, or Lgo. Both estimators 
are used under the assumption that all primary CR are 
either protons or iron nuclei. These extreme assumptions 
lead to a bias which then has to be corrected for. 
Using Ng and slope the energy of the primary cosmic- 
ray nucleus is reconstructed in two basic steps here 
(Lindner 1998b; Cortina et al. 1997; Plaga ct al. 1995): 
first slope (a measure of the distance to the shower max- 
imum) is combined with iVg to estimate the number of 
particles in the shower maximum which is proportional to 
the energy contained in the electromagnetic component of 
the EAS. In the second step a specific primary mass is as- 
sumed; with the assumption of primary proton (iron) we 
denote the methods as 1 (2). This allows to calculate the 
primary energy from the energy deposited in the electro- 
magnetic component. The following relation was used in 
our analysis: 

log{E[TeV]) = a ■ log{N,{max)) + b. (3) 

Here a,b were obtained from the discrete Monte Carlo 
data as 0.965,-2.545 (0.890,-2.010) for protons (iron). 
A^s(max) is the shower size at the maximum of shower 
development and is inferred from Ng as: 

iVs(max)/A^s — ao + aislope -f a2slope^ + a^slope^ (4) 

with a„ given as (0.57833,-85.146,6181.8,-71054) for ah 
primary nuclei. This procedure is valid because the shape 
of the shower development is only weakly dependent on 
the mass of the CR nucleus A, especially after the shower 
maximum (Lindner 1998a). Only the fraction of the total 



energy fed into the electromagnetic cascade depends on 
A for a given energy per nucleus. The comparison of the 
results assuming initially proton and iron primaries is a 
consistency check for the dependence of shower size at the 
maximum of shower development on the energy per nu- 
cleon. 

Alternatively the energy is reconstructed from the AIRO- 
BICC data alone (method 3 (4) with the assumption of 
proton (iron) primaries). Here it turns out that Lgo is a 
good estimator of the energy contained in the electromag- 
netic EAS cascades. Prom simulations the relation 

log{E[TcV]) a ■ /og(Lgo[photons/m2]) + b (5) 

was derived, where the coefficients a and b are given as 
0.958,-1.810 (0.840,-1.061) for primary protons (iron). 
slope (the parameter used to estimate the primary mass 
composition, see next section) is not involved in this en- 
ergy reconstruction. For ease of reference the four energy- 
reconstruction methods are summarised in Table 1. The 



Table 1. Summary of the energy reconstruction methods 1-4 
as discussed in the text. 



Method # 


parameter used 


primaries assumed to be 


1 


A^s and slope 


protons 


2 


Ns and slope 


iron nuclei 


3 


light density Lgo 


protons 


4 


light density Lgo 


iron nuclei 



agreement of analyses based on Ng and Lgg is a consis- 
tency test for the accurate description of the longitudinal 
shower development by the Monte Carlo simulation. 
Naturally (because the fraction of the primary energy de- 
posited in electromagnetic cascades depends on the energy 
per nucleon of the primary particle) the mean of the calcu- 
lated energy is only correct for the assumed particle type 
(Figure 2). The biases shown in Figure 2 have to be cor- 
rected for, to derive the real energy spectrum and CR mass 
composition from our measurements (see section 5.2, 5.3). 
In order to check that our final results do not depend on 
the assumed primary-particle mass, we shall always com- 
pare the results based on the four energy reconstruction 
methods below. 

The distribution of the reconstructed energy compared to 
the simulated energy is shown for examples in Figure 3. 
Note that the energy reconstruction from Lgg alone shows 
Gaussian distributions while the energy obtained from Ng 
and slope exhibits tails to high values which have to be 
taken into account properly in the analyses. Figure 4 shows 
the relative energy resolution achieved for different pri- 
mary particles and energy reconstruction methods 1 and 
3. If Ng is involved in the energy reconstruction the c;ncrgy 
resolution is limited by the experimental accuracy of the 
shower size determination at the detector level. Due to the 
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Fig. 2. The bias of the energy reconstruction as a function of 
primary energy for different primary masses. Shown is the ratio 
of the reconstructed energy with method 3, divided by the true 
energy (from the Monte-Carlo simulation). Very similar results 

are obtained from Ns and slope. The lines show fits used for 
convolution procedures to determine the final results (see text). 

smaller Ng of iron compared to proton induced EAS the 

accuracy of the energy reconstruction for iron showers is a 
little worse than for proton showers. The energy resolution 
obtained from Lgo is mainly determined by fluctuations in 
the shower development (being larger for proton than iron 
induced showers) and could not be decreased much by im- 
proving the detector. 

In all analyses below we bin the data in six equidis- 
tant energy intervals from logioi?reconstructed [TeV]=2.5 to 
logio^^reconstructed [TcV]=4.0 (scc e.g. Fig.9). Evcut Sam- 
ples defined to contain events in a certain rcconstnictcd- 
energy interval for the four energy-reconstruction methods 
then contain events with different true primary energies. It 
should always be kept in mind that these four samples are 
not independent because they are all based on the same 
total data sample. 

5.2. Chemical Composition 

The composition of CR is determined by analysing the 
EAS penetration depth (Xmax) distributions in inter- 
vals of the reconstructed primary energy. Information 
is contained in the differences of the mean Xmax val- 
ues for different primaries (protons penetrate about 100 
- 130g/cm^ deeper than iron in the energy range con- 
sidered here) and also in the the different fluctuations 
of the shower maxima position. Including experimental 



Fig. 3. Distribution of the ratio of reconstructed to MC gen- 
erated energy (300 TeV) for two primaries (left distribution: 

iron nuclei, right distribution: protons) and the two different 
energy reconstruction methods discussed in the text. Each dis- 
tribution is normalised to the same area. 



resolution we obtained RMS(Xi„ax,p) = 90 g /cm? and 
RMS(Xmax,i^e) = 5G.g/cm2 at IPeV). The RMS values 
of the depth distributions of Monte-Carlo events slightly 
decrease with rising energy, an effect that is partly due to 

an improving measurement of slope. 

We perform an analysis which uses both of these param- 
eters in one fitting procedure. As the error from such an 
analysis turns out to be already quite large, we do not per- 
form an analysis based on mean penetration depth alone. 
An analysis based mainly on the fluctuation of penetra- 
tion depths is discussed in section 7. 
The present data are not sensitive enough on the chemical 
composition to allow a analysis with several mass groups; 
therefore we restrict oiirsclves to a determination of the 
fraction of light nuclei (protons and helium) by fitting the 
expected to the measured X^ax distributions. To define 
the MC expectations for light nuclei, the generated distri- 
butions for primary protons and helium nuclei are added 
with weights of 40% and 60% (the ratio derived from direct 
measurements at energies around 100 TeV (Wiebel-Sooth 
et al. 1998)). The distribution of heavier nuclei is con- 
structed analogously by summing 65% oxygen and 35% 
iron induced EAS. Variations in this ratio at higher en- 
ergies are possible and are an additional potential source 
for systematic errors that is not further considered below. 
The spectrally weighted Monte-Carlo data are fitted to the 
measured penetration-depth distributions for each of the 
four energy-reconstruction methods. Because spectrally 
weighted Monte-Carlo data were available only for the en- 
ergy bins logio (-E^reconstructed) = 2.5 - 2.75 and 3.5 - 3.75 
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Fig. 4. The energy resolution obtained for different primary 
nuclei as a function of the generated MC energy. A^'s in brack- 
ets denotes an energy reconstruction that combines the mea- 
sured shower size at detector level and slope, L90 denotes the 
results obtained from the Chcrcnkov light density alone. The 
energy reconstruction was always performed assuming that the 
primaries are protons (referred to as methods 1 (stars in the 
Figure) and 3 (dots in the Figure) in the text; this fact is sym- 
bolised by the subscript "p" on the energy) . The "proton" (full 
symbols) "iron" (open symbols) after the colon indicate the 
primary for which the energy resolution was determined. The 
lines show fits used for convolution procedures to determine 
the final results (see text). 



Fig. 5. The differential 

shower-size spectrum and "light-density at 90 m core distance 
(L90)" spectrum. The values used for the construction of these 
spectra, were employed for the energy reconstruction. The full 
lines indicate the best fit in the range 5.3-6.8 and 4.5-6.1 for Ns 
and L90 respectively. Two power laws which meet in a single 
flux value ("the knee") were assumed. The best-fit power law 
indices are (-2.35/-2.92) and (-2.61/-3.13) (before/after) the 
knee at a position of logio(A*'s)/logio(I/9o) = 5.99/5.64 for the 
shower size/light density respectively. The energy scales wore 
derived under the assumption that the primaries are all pro- 
tons, resp. iron nuclei for shower size (upper scales) and light 
density (lower scales). 



the energy bins between 2.5 - 3.25 (3.25 - 4) were fitted 
with the former (latter) distribution. The MC events used 
in energy intervals other than the two for which the sim- 
ulations were done, were shifted in the mean penetration 
depth according to the elongation rate of the various el- 
ements. To avoid any systematic uncertainties related to 
imperfect parameterisations of the MC distributions and 
to take into account the statistical uncertainty of the sim- 
ulated event sample we directly fit the MC generated dis- 
tributions to the experimental data. 

Due to the primary dependent energy-reconstruction 
method the results for the "fraction of light nuclei" (ab- 
breviated "(p + a)/air' below) are biased. The results 
for these fits in the chosen energy bins are shown for 
method 3 in Fig. 9. The obtained (p -I- a)/all ratios are 
then corrected for the A dependent bias which is illus- 
trated in Fig. 2. The correction can be described as a sin- 
gle overall factor for the (p -|- a)/all ratio for each en- 
ergy bin - rather than a transformation of the penetration 
depth distribution - to a good approximation because of 



the independence of our energy reconstruction methods 
of Xmax as discussed in section 5.1. These correction fac- 
tors were derived from spectrally weighted Monte-Carlo 
data via determining the true (p -|- a)/all in the Monte 
Carlo that yields the fitted biased {p + a)/ all in the given 
reconstructed-energy bin. In this way the ratio of biased 
to true (p -I- a) /all at the true mean energy of the Monte- 
Carlo showers in the energy bin for a given energy recon- 
struction method is obtained. As an illustration the correc- 
tion factors for the case of energy reconstruction method 
3 are shown in Table 2. 

For the spectral weighing of the Monte-Carlo sampling a 
primary-spectrum as obtained from low-energy measure- 
ments (Wiebel- Sooth et al. 1998) with a power law index 
of «=— 2.67 and a "knee" at 3.4 PeV with a change in the 
power-law index to a=— 3.1 was assumed. An iterative 
repetition of this procedure with the energy spectrum as 
inferred below from the present data is possible. However, 
it was found that the contribution to the systematic error 
introduced by not performing the iterations is negligible 
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for the initial parameters chosen. 

Two Monte-Carlo samples were used for bias corrections 
in this work, the Monte-Carlo sample with events con- 
tinuously distributed in energy, mentioned in section 3, 
and a "toy Monte-Carlo sample" with unlimited statis- 
tics, which was created by randomly choosing all mea- 
sured parameters (like reconstructed energy, X^ax etc.) 
of a shower with a given true primary energy from one 
dimensional distributions inferred from the Monte Carlo 
sample with discrete energies. It was checked that the cor- 
rections obtained with these samples are very similar in 
energy regions where the continuously distributed Monte- 
Carlo data were available. 

Table 2. The fraction of light elements (uncorrected) and cor- 
rection factor for the A dependent bias with energy reconstruc- 
tion method 3. The uncorrected ratio has to be multiplied by 
this factor to yield the unbiased ratio. The energy intervals are 
specified as the logarithm to the base of ten in units of TeV. 



Reconstructed energy p -I- a/light correction factor 



2.5 - 2.75 0.502 0.741 

2.75 - 3. 0.493 0.763 

3. - 3.25 0.553 0.813 

3.25 - 3.5 0.406 0.800 

3.5 - 3.75 0.272 0.806 

3.75 - 1. ().()!) 0.820 



5.3. Energy spectra, elongation diagrams and penetration 
depth fluctuations 

Energy spectra obtained with 

the four energy-reconstruction methods were corrected for 
the A dependent bias by dividing the flux values in bins 
with true and reconstructed energy in the Monte-Carlo 
samples. The chemical composition as determined with 
the methods in the previous section is used. These factors 
were applied to the flux in each energy bin when going 
from reconstructed (Fig. 6) to true energy (Fig. 7). 
-'^max as a function of true energy is obtained if the mean 
-'^max is plotted at the mean true energy of the events in 
a given reconstructed-energy bin, as calculated with the 
measured chemical composition. This procedure leads to 
correct results as long as the elongation rate of different 
nuclei is identical; this is fulfilled to a good approximation 
for all hadron generators. 

The RMS of the shower penetration depth distributions 
were directly calculated from the distributions calculated 
with a given energy-reconstruction method, i.e. no pro- 
cedure to remove the bias was applied. These results 
were compared with RMS values from Monte Carlo data 
treated in the same way. 



5.4- Experimental Statistical and Systematic Uncertain- 
ties 

For the energy spectrum the statistical uncertainties cor- 
respond to the square root of the energy-bin contents N 
for the energy spectrum and the mean Xmax divided by 
\/N for the penetration depth. In all other cases statistical 
errors were obtained by changing the fit parameter from 
its best-fit value until the Xred increases by 1. In case of 
best fit Xred^ excess of 1.5 the best fit value of the fit 
parameter was increased until Xred doubled. 
Systematic uncertainties of the Monte-Carlo simulation 
of hadronic air-showers - estimated by using different 
hadronic Monte-Carlo generators - will be considered in 
a forthcoming paper. Here we concentrate on experimen- 
tal uncertainties related to the slope reconstruction. These 
are contributions from remaining uncertainties in the char- 
acteristics of the AIROBICC amplifier (3% uncertainty 
for slope) and non-perfect knowledge of the layer struc- 
ture and the light absorption of the atmosphere above 
the detector (4% and 2%). Models of the atmosphere 
have been carefully checked using the large statistics of 
photon induced air showers which were registered with 
the HEGRA imaging air Cherenkov telescopes in 1997 
(Konopclko ct al. 1999). Added in quadrature the system- 
atic uncertainty of slope amounts to 5%. The mean AT max 
for 300 TcV proton (iron) induced showers is then de- 
termined with an uncertainty in the absolute values of 
20 (13) g/cm^. The uncertainties for different primaries are 
strongly correlated. 

For the chemical composition, the energy spectrum and 
the variation of A'max with energy (elongation rate), the 
systematic error was evaluated by changing all slopes by 
5 % (the systematic error of this parameter) up or down. 
The whole analysis, including energy reconstruction, was 
then repeated and the deviation of the results thus ob- 
tained to the original ones was taken as the systematic er- 
ror (errors beyond the tick mark in Figs. 10,11 and shaded 
bands in Figs. 7,8). The shaded band in Fig. 11 is ob- 
tained by varying the best-fit composition within its total 
systematic and statistical error. In the case of the elonga- 
tion rate the systematic error was found to be dominated 
by the differences in the four energy-reconstruction meth- 
ods, this is dicussed in detail in the Results section 6.4. 
In case of the RMS of the penetration depths, the spec- 
tral fit parameters (knee position, power-law indices) and 
the elongation rate, the systematic error was estimated as 
the sample standard deviation of the best fit parameters 
obtained with the four energy-reconstruction methods. 



6. Results 

In this section the methods explained in section 5 are ap- 
plied to the data set discussed in section 4. 
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6.1. Ns and Lgo spectra 

Fig. 5 shows the Ng and Lqq spectra. These spectra display 
a relatively sharp knee at values consistent with a primary 
energy for the knee as determined below. 



6.2. Energy Spectra 

Figure 6 displays the integral energy spectra uncorrected 
for an A dependent bias obtained with the four recon- 
struction methods. The differences in absolute normalisa- 
tion and spectral slope originate from the different mass 
dependent biases. After the correction of the chemical bias 




Fig. 6. The integral energy spectra obtained with the four 
energy reconstruction methods denoted with different symbols. 
Pull square: Method 1, open square: method 2, full triangle: 
method 3, open triangle method 4. The biases of the energy 
reconstruction have not been corrected for. 



the integral spectra are similar (Figure?). This nontriv- 
ial fact is in favour of the internal consistency of data 
analysed here; a longitudinal shower development differ- 
ent from the one predicted by the Monte Carlo or errors 
in the calibration of Lgo and could have spoiled the 
agreement of spectra obtained with different energy re- 
constructions. The differential energy spectrum is shown 
in Figures. A steepening of the energy spectrum is vis- 
ible around an energy of 4PeV. There seems to be no 
"fine structure" in the energy spectrum around the knee 
in excess of 20 %. Apparent structure with smaller am- 
plitudes that appears in the spectrum reconstructed with 
a given enc;rgy-reconstruction method is not reproduced 
with other methods. This is expected due to the A de- 
pendent bias of our energy-reconstruction methods (see 
Fig. 2). Note that with these methods, a potential struc- 
ture in the energy spectrum consisting of different nuclei 
is smeared out. If two different power laws, smoothly con- 




Fig. 7. Integral cosmic-ray spectrum corrected for the A de- 
pendent bias. The shaded area denotes the systematic error. 
Symbols are the same as in the previous Figure. 



nected at the knee (corresponding to a "sharp" knee), are 
fitted to the differential spectra we obtain: 



"knee" 

j-l-4.66 



position 



of 



i;(Knee) =3.98Io.83(stat) ± 0.53(syst) PeV, 



- a spectral index of -2.72lQ 03(stat) ± 0.07(syst) below 
and -3.22tg;^^(stat) ± 0.08(syst) above the knee. 

The reduced values of the fits to the differential spec- 
trum (12 d.o.f.) were 6.75, 4.03, 3.53 and 1.47 with en- 
ergy reconstruction methods 1 to 4 respectively. Some of 
these values are much larger than one. It is then diffi- 
cult to specify a statistical error; we specify the statistical 
errors for method 4 that has a marginally acceptable re- 
duced value. The large values for the analysis with 
energy-reconstruction methods 1-3 can be interpreted as 
an argument in favour of a knee not absolutely sharp in 
energy. However, the fact that one of the fits is acceptable 
on the 90 % confidence level means that we cannot reject 
the hypothesis of a "sharp" knee (two power laws with no 
transition region) within our systematic errors. The large 
statistical error on the knee position further indicates that 
we cannot reject the hypothesis of a spectrum without a 
knee in the limited energy range of this analysis with high 
significance. 

The spectral index for the spectrum below the knee is 
consistent with direct measurements at lower energies 
(Wiebel-Sooth et al. 1998) and a recent Cherenkov-light 
based determination of the spectral index in the TeV range 
(Aharonian et al. 1999); there is therefore no evidence for 
any change in spectral index from the TeV range right up 
to the knee. 
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Fig. 9. The fit of MC expectations for light and heavier nuclei to the measured shower maximum depth distribution in the 
analysed reconstructed-cnorgy intervals for method 3. The numbers in the upper loft corner arc the logarithms to the base of 
ten of the energy-bin boundaries in TeV. The full dots mark the experimental data with statistical errors, the crosses with error 
bars are the fitted MC distribution (where the errors correspond to the MC statistics). The two components fitted to the data 
are shown as dark shaded (large penetrations depths, light nuclei) and light shaded (heavy nuclei) histograms. Details of the 
procedure are described in the text. 



6.3. Composition of CR 

The fraction of light nuclei as a function of recon- 
structed energy - obtained from the fits to the measured 
penetration-depth distributions (see section 5.2)- is pre- 
sented in Figure 10. At energies below the knee the com- 
position is mixed and consistent with direct measure- 
ments around 100 TeV, namely {p + a)/sl\ = 0.54 ± 0.08 
(Watson 1997). The data points seem to indicate a gradual 
enrichment of heavy elements above about 1 PeV though 
the error bars are large (remember that there are only 
six independent data points). We will discuss in section 



7 how reliable the qualitative conclusion of a gradual en- 
richment in heavy elements is within our systematic errors. 
The data rule out a predominantly light composition at all 
energies and does not give evidence for a drastic change 
of composition at the knee. 



6.4- Elongation rate 

Figure 11 shows the corrected mean shower maximum 
depth as a function of energy. A least-squares fit to the 
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g/cm^ from 300 TeV up to 10^ TeV (dotted fit line in the 
summary diagram 10 in (Watson 1997)) is an intriguing 
hypothesis which is not in contradiction with our data. 

6.5. Fluctuation of shower penetration depth 

The RMS of the penetration depth distributions - calcu- 
lated in reconstructed-energy bins, i.e. biased in favour of 
the light component of CR especially at low energies - is 
shown in Table 3. It does not show any obvious trend to- 
wards a heavy composition. Therefore the fact that the 
composition at the highest energy seems to be heavy with 
all energy reconstruction methods (Fig. 10) is mainly de- 
termined by the fact that the Xmax in the highest energy 
bin lies about 20 g/cm^ below a constant elongation rate. 



Fig. 8. The differential CR energy spectrum obtained with 
four energy reconstruction methods. Symbols are the same as 
in Fig. 6. The light shaded region represents systematical un- 
certainty. The "star" with vertical error bars shows the uncer- 
tainty of the HEGRA data originating from the 10% systematic 
uncertainty of the absolute energy scale that we estimate from 
the uncertainty in the determination of the absolute Ns scale. 
The full line is the best fit as described in the text. 

Xaia.x values as a function of energy, using only the statis- 
tical errors, 



X„,ax = ERx logio(E) + ERB 



(6) 



yields an elongation rate ER=78.3 ± 1.0 (stat) ± 6.2 
(syst) g/cm^ and mean depth parameter ERB=243.1 ± 
2.6 (stat) ± 15.7 (syst) g/cm^. The specified mean values 
and statistical errors are the mean of fit values with the 
four energy-reconstruction methods. The systematic error 
is estimated as the standard deviation of the mean val- 
ues inferred with the four energy-reconstruction methods. 
The systematic error introduced by the systematic uncer- 
tainty in slope is smaller (about 3 and 14 g/cm^ for ER 
and ERB respectively). The reduced values of the fit 
to relation (6) (4 d.o.f.) are very large (6.6,9.2,17.2,23.5) 
for energy-reconstruction methods 1-4, i.e the systematic 
errors dominate over the rather small statistical errors for 
the mean Xmax- Therefore the specified estimates of the 
statistical errors obtained with the procedure explained 
in section 5.4 have to be treated with caution. The data 
point at the highest energy lies about 20 g/cm^ higher in 
the atmosphere than expected for a constant elongation 
rate. 

These results are not in contradiction with previous 
measurements in this energy range (Wdowczyk 1994; 
Turver 1992). This elongation rate, and also the abso- 
lute X„iax, is consistent with data at higher energies, ob- 
tained mainly by the Yakutsk and Fly's Eye collabora- 
tion (Watson 1997). A constant elongation rate of ^ 73 
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Fig. 10. The corrected fraction of light nuclei determined with 
four different energy reconstruction methods (see text for de- 
tails). Each data point is plotted at the true mean energy of 
the events used to infer the mean Xmax- The error bars are 
statistical and systematic error added in quadrature (up to 
the tick- mark: only statistical error) . The statistical errors are 
correlated due to the use of an identical Monte-Carlo sample 
in the first and last three energy intervals. The shaded band 
shows the allowed region between a polynomial fit to the upper 
and lower ends of the error bars. 



7. Further studies of systematic uncertainties; 
analysis methods independent of absolute X^ax 

The trend for an enrichment in heavy elements above the 
knee - which Fig. 10 suggest - is not significant within our 
errors. In this section we elucidate this fact further, and 
explore what it would take to detect significantly a modest 
trend for an enrichment in heavy elements - as expected 
e.g. in a diffusion model of the knee (see Introduction) - 
with the present techniques. The agreement between the 
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Table 3. The RMS of the penetration depth distributions [g/cm'^] as a function of reconstructed energy (given in the same 
units as in Table 1) in the data and spectral Monte-Carlo sample. Given is the value inferred for the energy bins as defined 
with energy-reconstruction method 3, i.e. the specified values contain an A dependent bias. The first error is statistical and the 
second systematic (due to the systematic error in slope). For the numbers from Monte-Carlo simulations only a statistical error 
is given. "Mixed composition" represents the expectation for our best-fit chemical composition. Based on numerical experiments 
the statistical error was taken as the inferred value divided by V number of events in bin rather than half of that, as would be 
correct for Gaussians, due to the non-Gaussian tails of the distribution. The comparison between experimental data and Monte 
Carlo simulations shows no trend towards a heavy composition at the higher energies. 



Rec. energy Data MC:mixcd comp. MC:p MC:Fe 

2.5 - 2.75 (ca. 0.32 - 0.56 PeV) 84 ± 0.6 ± 1.8 77 ± 2 83 ± 4 50 ± 3 

2.75 - 3. (ca. 0.56 - 1 PeV) 80 ± 0.9 ± 1.4 

3. - 3.25 (ca. 1 - 1.79 PeV) 80 ± 2 ± 1.0 

3.25 - 3.5 (ca. 1.79 -3.16 PeV) 73 ± 2 ± 1.1 

3.5 - 3.75 (ca. 3.16 - 5.62 PeV) 67 ± 3 ± 1.5 45 ± 5 73 ± 17 48 ± 11 

3.75 - 4. (ca. 5.62 - 10 PeV) 67 ± 7 ± 0.9 
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Fig. 11. The mean shower maximum depth as a function of 
energy using QGSJET simulations to model the EAS develop- 
ment. To obtain an unbiased elongation plot each data point is 
plotted at the true mean energy of the events used to infer the 
mean Xmax- The shaded region indicates the region expected 
from our best fit composition within its total error. Errors are 
statistical and systematical errors added in quadrature, up to 
the tick mark only statistical. Up to the highest energies the 
systematical error dominates. 



Xmax distribution shape at low energies - predicted as- 
suming a composition at low energies which is not very 
different from the one obtained by direct experiments - 
and the data (see Fig. 9) is satisfactory. This is an argu- 
ment in favour of a correct MC simulation. 
To explore the effect of our systematic error in slope (sec- 
tion 5.4) Fig. 12 shows the results derived with an initial 
assumption of slope changed by 5 % from its preferred 
value for energy-reconstruction method 3 and 4. It be- 
comes clear that not only the mean but even the over- 
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Fig. 12. The inferred chemical composition using en- 
ergy-reconstruction method 3 (upper panel) and 4 (lower 
panel). The dots are the values as discussed before, the squares 
(triangles) are the results obtained when the slope is increased 
(decreased) by 5 % (the systematic error on this variable) . The 
general "trend" (composition gets heavier/lighter) may change 
within this systematics. 



all apparent "trend" may change, e.g. for method 3 with 
slope increased by 5 % the composition appears to become 
lighter from the knee up to the penultimate bin. It should 
be stressed that none of the discussed "trends" is signifi- 
cant within our statistical errors. 

The deviation of the penetration depth at the highest en- 
ergy point from a constant elongation rate discussed in 
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subsection 6.5 is of the order of disagreements between dif- 
ferent Monte Carlo codes at this energy (Heck et al 1998). 
Therefore the possibility remains open that this deviation 
is due to a change in cascade characteristics not reliably 
modelled by the Monte Carlo, rather than to an enrich- 
ment in heavy elements. As it is also quite similar in size 
to our systematic error in slope, an origin in the HEGRA 
experiment for this deviation is also difficult to rule out. 




Fig. 13. The reduced ^ values (z-axis) of a description of 
the measured penetration- depth distribution with the spec- 
tral Monte Carlo data as a function of "fraction of light nu- 
clei" (y-axis) and overall shift in depth (x-axis). Only ac- 
ceptable values smaller than 1.5 are displayed. The en- 
ergy was reconstructed assuming protons and using the the 
Cherenkov light density (method 3). The left panel is for 
the first energy bin (logioE[TeV]=2.5-2.75, 19000 data events, 
1460 Monte-Carlo events) the right panel for the penultimate 
one (logioE[TeV]=3.5-3.75, 369 data events, 98 Monte-Carlo 
events) . In the high-energy bin practically all chemical compo- 
sitions are allowed for certain "shifts". 

Previous conference publications (Plaga et al. 1995; 
Cortina et al. 1997a; Cortina et al. 1998) are superseded 
by the present results - differences are mainly due to a 
more sophisticated amplifier calibration and the simpler 
energy reconstruction for the present analysis. In the last 
two of these publications, we tried to lessen the depen- 
dence of our composition result on the correct absolute 
-''^max values by allowing a free "shift parameter" ; the fit 
to the Xniax was then performed with two free parame- 
ters: the ratio of light to all nuclei and an overall shift in 
penetration depth of all MC distributions. In this way the 
result is mainly determined by the shape of the Xmax dis- 
tributions (in first order its width, i.e. RMS value). This 
width depends only weakly on a systematic uncertainty in 
the determination of slope, relative to the expected differ- 
ence of a purely light or heavy composition. 
Fig. 13 displays the result of such an attempt in a two di- 
mensional plot showing the reduced for various "frac- 
tion of light nuclei" - "shift parameter" combinations for 



the data lowest and highest in energy. The shift is varied in 
an interval ± 30 g/cm^, estimated from the likely system- 
atic uncertainty of our detector and the Monte Carlo code. 
While in the low energy bin small (p-f-a) /all ratios lead to 
unsatisfactory values for all shift values, in the highest 
energy bin - well above the knee - practically all fractions 
give acceptable values for appropriate shifts as seen 
in Fig. 13. The reason for this behaviour is that - given 
the small number of events in the high-energy bin - the 
Xniax distribution can be fitted both with the relatively 
broad predominantly light composition shifted to larger 
depths in the atmosphere and a mixed heavy/ light com- 
position (where the difference in mean penetration depth 
of the heavy and light component contributes to the to- 
tal width) shifted to small penetration depth. We have to 
conclude that it is not reliably possible to determine the 
composition based mainly on the width of the X^ax distri- 
bution. Wc found in numerical experiments that with this 
method, and assuming a Monte-Carlo simulation describ- 
ing the experimental data well, together with a statistics 
increased by about a factor of 100 (which is difficult but 
not impossible to reach in future experiments) it will just 
be possible to reach the desired precision of 10 % men- 
tioned in the introduction on a 2a level beyond the knee. 



8. Conclusion 

The results of this paper demonstrate that our data seem 
to be well described by the QGSJET hadron generator 
within the CORSIKA program: The energy spectra de- 
rived with different air-shower observables and assump- 
tions about the chemical composition of CR differ only 
within the estimated systematic errors. This is an argu- 
ment in favour of the contention that the "knee" is a fea- 
ture in the primary flux of cosmic rays, rather than some 
new effect in the interaction of cosmic rays at very high 
energies. Such a new effect would probably lead to incon- 
sistencies in the data analysed by a generator that does 
not take them into account. 

Up to the knee we find a mixed composition and an 
energy-spectrum power law index consistent with the re- 
sults of direct experiments at energies around 100 TeV. 
The data favour a gradual enrichment in heavy elements 
at energies above the knee but the systematic errors of our 
experiment do not allow to rule out a constant composi- 
tion. An abrupt change of composition or a substantial 
enrichment in light elements at the knee is ruled out. We 
find an elongation rate of about 78 g/cm^ and a mean pen- 
etration depth consistent with other experimental data on 
CR with energies higher than studied here (Watson 1997). 
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